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Broadband  Spectral  Estimation  Using  The  Multiple  Window  Method 
A  Comparison  With  Classical  Techniques 

Jean-Marie  Tran 

Marine  Physical  Laboratory 
Scripps  Institution  Of  Oceanography 
San  Diego,  CA  92152 


ABSTRACT 


The  Multiple  Window  recently  introduced  by  D.J.  Thomson  is  an  alternative 
approach  to  spectral  estimation.  Using  an  orthogonal  decomposition  of  the  Discrete 
Fourier  Transform  of  the  data  on  the  prolate  spheroidal  wave  sequences,  the  power  spec¬ 
trum  is  estimated  by  a  weighted  average  of  raw  spectra  obtained  from  the  Fourier 
transforms  of  the  data  windowed  by  the  prolate  sequences.  This  method  has  the  reputa¬ 
tion  of  working  well  on  short  data  sequences,  offering  good  variance  control  and  a  fair 
resolution.  Because  the  Multiple  Window  technique  combines  Fourier  transforms  with 
data  adaptive  weighting,  it  is  compared,  in  a  Monte  Carlo  simulation,  to  more  conven¬ 
tional  spectral  estimation  methods  :  the  classic  averaged  Fourier  transforms  technique 
and  Capon’s  Maximum  Likelihood  method.  The  underlying  spectral  density  to  be 
estimated  is  of  an  ARMA  process  simulating  underwater  acoustic  ambient  noise.  It  is 
characterized  by  a  moderate  dynamic  range  of  30  dB  and  can  be  regarded  as  a  generic 
spectrum  of  many  applications.  The  three  methods  considered  here  are  non  parametric, 
so  that  no  estimator  has  an  unfair  advantage.  The  normalized  bias,  random  and  root  mean 
square  errors,  as  well  as  the  statistical  distributions,  based  on  512  realizations  of  the  spec¬ 
tral  density,  quantify  the  performance  of  the  three  different  techniques.  / 


1.  Introduction 


The  multiple  Window  method  introduced  by  Thomson  [1,2,3]  is  a  new 
spectral  estimation  method  that  has  the  reputation  of  working  well  on  short  data 
sequence,  offering  good  variance  control,  bias  protection  and  fair  resolution. 

It  is  of  interest  to  compare  in  a  Monte  Carlo  simulation  the  resolution  capa¬ 
bilities  and  the  statistical  characteristics  of  the  Multiple  Window  method  with 
classical  methods  such  as  the  conventional  FFT  spectral  estimation  method  and 
the  Maximum  Likelihood  method. 

This  simulation  focuses  on  broadband  spectral  estimation.  The  reference 
spectrum,  estimated  by  the  various  methods,  is  characteristic  of  an  underwater 
acoustic  ambient  noise  spectrum  between  DC  and  250  Hz.  The  main  features  of 
ambient  noise  spectra  are  a  moderate  20  to  30  dB  dynamic  range,  a  red  part  at 
low  frequency  and  an  almost  white  part  at  high  frequency  [4],  This  reference 
spectrum  can  be  regarded  as  a  generic  spectrum  characteristic  of  many  applica¬ 
tions. 

After  an  outline  of  the  analysis,  some  background  material  relevant  to  the 
Multiple  Window  method  and  the  Maximum  Likelihood  method  are  presented. 
Then  the  ARMA  process  creating  the  data  that  simulate  the  ambient  noise  is 
described.  Finally,  results  of  the  simulation  are  given  and  discussed. 


2.  Analysis  Description 


The  sampling  frequency  is  assumed  to  be  500  Hz  so  that  the  power  spec¬ 
trum  lies  between  DC  and  250  Hz.  All  the  results  presented  on  the  plots  are 
given  in  terms  of  the  fractional  sampling  frequency,  between  0  and  0.5  of  the 
Nyquist  frequency. 

Each  spectral  estimate  is  characterized  by  the  following  plots  : 

(1)  the  overlaid  realizations  of  the  estimated  power  spectrum 

(2)  the  mean  and  envelope  of  the  realizations  of  the  power  spectrum 

(3)  the  normalized  random  error 

(4)  the  normalized  bias  error 

(5)  the  normalized  root  mean  square  error 

(6)  the  distributions  of  the  estimate  at  the  normalized  frequencies  :  10/512, 
50/512,  200/512  (or  9.76  Hz,  48.8  Hz  and  195  Hz),  surimposed  with  fitted 
Gaussian  and  chi-square  distributions  converted  in  number  of  observations. 

The  error  estimates  allow  an  overall  rating  of  the  different  methods.  The  nor¬ 
malized  errors  are  defined  as  in  [5],  The  normalized  random  error  er  is 


^VAR(<t>) 


(1.1) 


where  VAR  is  the  variance  operator,  is  the  power  spectral  density  estimate  and 
4>  is  the  true  underlying  spectrum.  The  normalized  bias  error  eb  is 


E(<t>) 


where  E  is  the  mean.  The  root  mean  square  (rms)  error  e^,  is 
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The  three  selected  frequencies  correspond  to  the  low  frequency  (highly 
colored),  intermediate  frequency  (pink)  and  high  frequency  (almost  white) 
behaviors.  The  distributions  give  additional  information  on  the  statistics  of  the 
estimates  and  are  therefore  of  interest.  The  Gaussian  distribution  is  plotted  for 
reference  and  is  entirely  determined  by  the  sample  mean  and  variance  of  the 
spectral  estimates  at  a  selected  frequency  bin.  The  chi-square  distribution  is 
characterized  by  the  mean  of  the  estimates  and  the  number  of  degrees  of  free¬ 
dom.  The  analytical  expression  for  the  chi-square  with  mean  p.  and  order  n  is 
given  by  [5] 


’(*)  =  22rt|)  ^ 


where  r  is  the  Gamma  function. 

Otherwise  stated,  the  selection  of  the  fitted  chi-square  distribution  is  based, 
in  each  case,  on  the  chi-square  goodness  of  fit  test  [5]  using  10  constant  intervals. 
The  discrepency  between  the  sample  histogram  and  the  fitted  chi-square  is  meas¬ 
ured  by 

to  [/.  -/=•.]  2 

X2=S  -  -  ~  (1-5) 

i»l  r‘ 

where  /,  is  the  observed  number  of  observations  and  F,  the  number  of  observa¬ 
tions  corresponding  to  the  fitted  distribution  in  the  i*  interval. 

The  acceptance  region  of  the  fit  between  the  sample  histogram  and  the  chi- 
square  distribution  at  the  a  =  0.05  level  of  significance  is 

*2  5  X2*,o.o5  (1.6) 

The  number  of  degrees  of  freedom  for  which  %2  is  calculated  is  8  since  there  are 
9  independent  frequencies  /,  and  the  mean  of  the  fitted  distribution  is  estimated 
by  the  sample  mean.  The  chi-square  goodness  of  fit  therefore  produces  a  rough 
estimate  of  the  number  of  degrees  of  freedom  of  the  estimates  which  is  a  simple 
characterization  of  the  distribution,  assuming  it  follows  a  chi-square  law. 

3.  Outline  Of  The  Multiple  Window  Method 

The  multiple  window  (MW)  method  is  described  at  length  in  [1,2,3].  It  is 
not  the  purpose  of  this  section  to  give  a  complete  mathematical  derivation  of  the 
method  which  starts  with  the  Cramer  spectral  representation  of  random 
processes.  Rather  the  goal  of  this  paragraph  is  to  give,  in  a  short  summary,  back¬ 
ground  information  on  the  MW  method. 

The  MW  method  is  a  non  parametric  spectral  estimation  method  based  on 
multiple  windowed  Fast  Fourier  Transforms.  The  computation  of  the  multiple 
windows,  also  called  the  Slepian  windows,  is  a  preliminary  step.  The  method 
uses  a  subset  S  of  orthonormal  windows  that  are  the  eigenfunctions  of 

(21) 

where  the  bounds  of  the  integral  can  be  either  -Vi ,  Vi  or  -W,W  (W  is  the  half 
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bandwidth  of  the  window,  a  frequency  parameter  that  characterizes  the  MW  esti¬ 
mates,  given  N  data  points),  y(f)  is  the  Discrete  Fourier  Transform  of  the  data 
sequence  of  length  N,  dZ(y)  is  the  unknown,  related  to  the  power  spectral  density 
S(f)  by 

E[  \dZ(f)\2]  =  S(f)df  (2.2) 

The  driving  idea  lies  in  the  expansion  of  the  solution  dZ(y)  in  term  of  the 
eigenfunctions.  The  subset  S  includes  the  2NW  windows  that  have  the  most  con¬ 
centrated  energy  in  the  bandwidth  2W  (it  is  understood  that  2NW  is  the  integer 
part  of  2NW),  the  energy  content  of  a  window  corresponds  to  its  eigenvalue,  Xt, 
associated  to  the  integral  equation  (2.1).  Bias,  which  can  be  thought  of  as  noise, 
increases  as  one  keeps  higher  and  higher  order  windows  with  lower  and  lower 
energy  and  can  be  evaluated  by  its  upper  bound  given  by  the  Schwartz  inequality 

Bk  Or)  =  o2(l-X*)  (2.3) 

where  cr2  is  the  power  of  the  data  sequence. 

The  2NW  Discrete  Fourier  Transforms  of  the  data,  windowed  by  each  one 
of  the  2NW  windows,  are  then  computed  ;  the  inferred  spectra,  Sk(f),  are  called 
eigenspectra.  The  Slepian  windows  give  limited  statistical  consistency  to  the 
eigenspectra  ;  windowing  is  interpreted  in  the  frequency  domain  as  an  averaging 
over  the  bandwidth  of  the  window,  that  is  2W. 

To  get  better  variance  control,  the  eigenspectra  are  finally  averaged  using 
adaptive  weights  that  minimize  the  mean  squared  error  of  the  spectral  estimate. 
The  final  estimate  §(f)  is  solution  of 

2NW-1 

X  — - f  .  L  .  - —  =  0  (2.4) 
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Because  of  the  orthonormality  of  the  windows,  degrees  of  freedom  are  indeed 
added,  thus  producing  a  consistent  estimator. 

The  method  is  easily  extended  to  mixed  spectra  using  statistical  tests 
OFisher’s  statistic  test)  to  detect  line  components  in  the  spectrum.  The  spectrum 
is  corrected  or  reshaped  once  a  line  component  is  detected.  The  power  in  the 
line  component  is  subtracted  from  the  spectrum  over  a  2W  wide  frequency  zone 
around  the  sinusoid  by  using  a  local  expansion  of  the  spectrum  on  the  Slepian 
windows. 

The  estimates  of  the  power  spectrum  discussed  hereafter  corresponds  to 
(2.4)  which  is  implemented  in  the  Fortran  programs  written  by  Allan  D.  Chave. 
These  programs  are  such  that  they  yield  correct  power  for  white  noise  [2].  This 
paper  only  deals  with  continuous  spectra,  although  one  of  the  strength  of  the 
MW  method  is  its  ability  to  handle  mixed  spectra. 

As  pointed  out  in  [1,2,3]  a  good  compromise  between  variance,  bias  and 
resolution  is  achieved  for  a  time  bandwidth  product  NW  of  4.  The  MW  method 
has  the  reputation  of  working  well  on  short  data  sequence,  offering  good  vari¬ 
ance  control,  good  bias  protection  and  fair  resolution.  To  quantify  the  MW 
characteristics,  the  selected  data  length  for  the  numerical  simulations  is  N=128. 


4.  Outline  Of  The  Maximum  Likelihood  Method 

As  mentioned  in  §3,  the  Multiple  Window  method  combines  with  classical 
Fast  Fourier  Transforms  a  data  adaptive  technique  to  weight  the  eigenspectra  and 
produce  the  final  spectral  estimate.  It  is  thus  of  interest  to  compare  the  MW 
results  with  those  of  a  classic  data  adaptive  estimation  method  such  as  the  Max¬ 
imum  Likelihood  method  also  called  the  Minimum  Variance  method  [6,7,8]. 

The  Maximum  Likelihood  (MLM)  method  designs  an  optimal  finite 
impulse  response  filter  that  minimizes  the  output  power  or  variance  and  passes 
the  center  frequency  undistorted.  The  power  spectrum  is  given  by 


where  R is  the  inverse  of  the  autocorrelation  matrix  of  order  p ,  H  the  Hermitian 
operator  and  e  is  the  complex  sinusoid  vector 

«</)=[  i  e/2*/r  . .  •  '  (3.2) 

t  represents  the  tranpose  operation,  /  the  frequency  and  T  the  sampling  period. 

Several  computational  methods  can  be  used  :  a  direct  inversion  of  the  autocorre¬ 
lation  matrix  (considering  the  possibility  of  the  matrix  being  singular  as  decribed 
in  [6]),  or  an  indirect  determination  of  the  power  spectrum  using  a  relationship 
between  the  Maximum  Likelihood  estimate  and  the  AR  spectral  estimate  [7,8]. 
The  method  used  in  this  simulation  computes  an  estimate  of  the  autocorrelation 
matrix  and  inverts  it.  A  check  for  singular  matrix  is  performed.  The  singular 
matrices  can  be  stabilized  by  adding  a  small  number  to  the  main  diagonal. 

The  power  spectrum  is  then  computed  by  Fast  Fourier  Transforms,  on  the 
basis  of  equation  (3.1)  rewritten  as 

eHR~'e  =  (eH  (R-1)')’ e‘  (3.3) 

where  *  is  the  complex  conjugate.The  sequence  of  operations  is 

(1)  the  inversion  of  the  autocorrelation  matrix  R  (a  Toeplitz  matrix) 

(2)  the  transposition  of  R 

(3)  the  FFT  of  all  the  columns  of  the  matrix  R~l  giving  a  (NFFT,p)  matrix 

(4)  the  complex  conjugaison  of  the  (NFFT, p)  matrix  of  step  3 

(5)  the  FFT  on  all  the  lines  of  the  (NFFT,p)  matrix  of  step  4,  giving  a  (NFFT,  NFFT) 
matrix 

The  power  spectrum  is  recovered  by  taking  the  inverse  of  the  real  part  of 
each  element  of  the  main  diagonal  of  the  ( NFFT  ,  NFFT)  matrix  of  step  5.  It  is  then 
normalized  by  the  sampling  time,  multiplied  by  2  to  to  yield  the  proper  one-sided 
spectrum  levels  and  finally  multiplied  by  p,  the  order  of  the  autocorrelation 
matrix  R ,  to  yield  proper  power  for  a  white  noise  [8]. 


i 


5.  Ambient  Noise  Data 

An  ARMA  process  of  order  (l,n)  is  used  to  model  the  ambient  noise  and 
create  the  simulation  data.  Typical  ambient  noise,  as  given  in  [4],  can  be  approx¬ 
imated  in  a  straighforward  manner  by  a  transfer  function  of  the  form 


1+1  atz' 


-  (4.1) 

1+1  bki~* 

By  visual  selection  of  the  poles  and  zeros  of  this  transfer  function  on  the  z  plane, 
a  theoretical  frequency  response  close  in  shape  to  an  ambient  noise  spectrum  can 
be  created. 

Using  an  iterative  process,  five  poles  were  selected  at 

z  =0.99 

2  =  0.7886876±  i  0.452389 


and  four  zeros  at 


2  =  0.899868±  i  0.015707 


2  =0.895456 ±i  0.090326 
2  =  0.8966357+ 1  0.080891 


Poles  and  zeros  are  respectively  plotted  on  Figure  5.1. 

The  corresponding  ARMA  process  of  order  (4,5)  has  the  following  transfer 
function 


«(*)  = 


1-3.5836  2~ 1  +  4,8306  z~2  -  2.90272  2~3  +  0,6561  2^ 

1  -  4.3671  2*1  +  7.8188  z~2  -  7.1962  z~*  +  3.4074  2^  -  0.66291  z~ 5 


The  corresponding  frequency  response  is  plotted  on  Figure  5.2. 

The  driving  sequence  is  a  centered  white  Gaussian  noise  with  unit  variance. 
The  ARMA  process  is  readily  implemented  as  finite  impulse  response  and 
infinite  impulse  response  filters. 

1024  points  of  the  data  sequence  are  plotted  on  Figure  5.3.  This  data  sample 
corresponds  to  2.05  seconds  with  a  sampling  frequency  of  500  Hz.  The  theoreti¬ 
cal  one  sided  "ambient  noise”  spectral  density  is  plotted  between  DC  and  250  Hz 
on  Figure  5.4. 
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Figure  5.3  Simulated  noise  data  sequence  (1024  points) 


Figure  5.4  Theoretical  "ambient"  noise  power  spectral  density 
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6.  Averaged  Periodogratn  Estimate 

The  reference  periodogram  estimates  are  incoherent  averages  of  eight  FFTs 
with  a  Hanning  window,  the  transforms  are  performed  on  128  data  points  with  an 
overlap  of  50  %.  The  length  of  the  FFT  is  512  (i.e.  128  data  points  zero-padded 
out  to  the  FFT  length  of  512).  Each  periodogram  estimate  thus  is  based  on  576 
data  points. 

A  fairly  small  data  length  is  used  in  the  transform  because  of  the  MW 
method  capability  to  work  on  short  data  sequences.  The  number  of  averaged 
FFTs  corresponds  to  the  numerical  workload  of  the  MW  method  when  the  time 
bandwidth  product  NW  is  4  and  the  number  of  windows  2NW  is  8. 

512  different  spectra  were  computed  on  the  simulated  ambient  noise  time 
series  and  are  overlaid  on  Figure  6.1.  The  envelope  and  the  mean  of  these  512 
estimates  are  plotted  on  Figure  6.2.  The  mean  gives  an  idea  of  the  estimate  bias 
with  respect  to  the  theoretical  ambient  noise  spectrum  (Figure  5.4).  The  envelope 
gives  an  idea  of  the  scatter  or  variance  of  the  FFT  based  estimates.  The  ambient 
noise  spectrum  is  well  resolved  by  this  classic  method. 

The  random  error  er  is  plotted  on  Figure  6.3  and  varies  between  0.35  and 
0.40.  The  spectral  estimate  is  the  average  of  nd  raw  spectra.  It  approximately 
follows  a  chi-square  distribution  with  2nd  degrees  of  freedom.  The  random  error 
is  evaluated  by  the  ratio  of  the  standard  deviation  of  the  estimate  to  its  mean  [5] 
giving  1  /  VnJ  ;  in  the  present  case  the  random  error  should  theoretically  be  0.35. 

The  bias  error  zb  plotted  on  Figure  6.4  is  low,  below  0.2,  except  at  low  fre¬ 
quencies  where  it  peaks  up  to  2.75.  The  rms  error  e**,  is  dominated  by  the  ran¬ 
dom  error.  It  is  plotted  on  Figure  2.5  and  is  bounded  by  0.35  and  0.4  except  at 
low  frequency  where  it  goes  up  to  3.46. 

As  mentioned  before,  the  distribution  of  the  estimates  is  theoretically  a  chi- 
square  distribution  with  16  degrees  of  freedom.  Fitted  Gaussian  and  fitted  chi- 
square  distributions  are  compared  to  the  histograms  of  the  spectral  estimates  at 
the  selected  frequencies  (§2).  They  are  respectively  plotted  on  Figures  6.6,  6.7 
and  6.8.  There  is  indeed  good  agreement  in  term  of  the  number  of  degrees  of 
freedom  of  the  estimate  since  they  are  respectively  16,  16  and  14  at  the  normal¬ 
ized  frequencies  10/512,  50/512  and  200/512. 
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Figure  6.1 


512  overlaid  conventional  FFT  spectral  densities  (8  incoherent 
averaged  FFTs  on  128  points  with  a  Hanning  window) 


Figure  6.2 


Mean  and  envelope  of  the  512  spectral  density  realizations 


Figure  6.3  Normalized  random  error  er 


Figure  6.4 


Normalized  bias  error  eb 


Figure  6.5 


Normalized  rms  error  e„ 


Figure  6.6 


Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  x2\e  distributions 


Figure  6.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x2is  distributions 


Figure  6.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2u  distributions 


7.  Multiple  Window  Estimate 

7.1.  Selection  Of  The  Time  Bandwidth  Product  NW 


As  stated  earlier,  the  optimal  time  bandwith  product  NW  is  4.  If  the  esti¬ 
mate  is  based  on  128  points,  the  equivalent  MW  filter  normalized  bandwidth  2W 
is  0.0625.  The  main  features  of  the  simulated  noise  is  a  -25  dB  dip  at  the  normal¬ 
ized  frequency  0.02  (or  8  Hz)  which  is  followed,  15  dB  higher,  by  a  moderate  Q 
resonance  (shipping  bump)  just  before  the  normalized  frequency  0.01  (or  around 
50  Hz).  One  can  already  predict  that  these  features  will  be  smeared  by  the 
equivalent  filtering  of  the  multiple  windows,  thus  leading  to  a  significant  bias  of 
the  estimate  at  low  frequency,  if  such  a  time  bandwidth  product  is  selected. 

To  illustrate  this  point,  64  spectra  were  computed  on  the  first  128  points  of 
each  periodogram  576  point  long  sequence.  NW  is  4  and  the  number  of  windows 
is  8.  No  reshaping  was  performed  since  there  is  no  line  component  in  the  data. 
The  smearing  effect  is  clearly  seen  on  Figure  7.1.1  and  Figure  7.1.2  where  the 
overlaid  spectra  and  the  mean  and  envelope  of  the  64  estimates  are  plotted. 

The  very  red  spectral  zones  are  smeared  near  DC  where  the  spectrum  is  flat 
up  to  a  frequency  that  corresponds  to  the  half  averaging  filter  bandwidth  W.  The 
resonance  is  broader,  with  a  lower  Q  than  the  theoretical  one. 

When  the  time  bandwidth  product  is  4,  the  averaging  equivalent  filter 
bandwidth  is  too  large  to  be  able  to  resolve  the  structure  of  the  simulated 
ambient  noise.  In  order  to  see  more  details  in  the  spectrum  estimate,  one  has  to 
decrease  the  time  bandwidth  product  but  the  drawback  is  more  variance  in  the 
estimate  [1,2,3]. 

The  analysis,  described  in  §2  and  carried  out  in  §5,  is  repeated  for  two  time 
bandwidth  products  NW  =  2  and  NW  =  3.  Then,  to  compare  the  classical  FFT 
based  method  and  the  MW  method  on  the  basis  of  the  same  information,  MW 
estimates  are  computed  on  the  same  576  point  sequences  as  in  §6,  with  a  time 
bandwidth  product  NW  of  4  and  8  windows. 


64  overlaid  spectral  densities  computed  by  the  MW  method 
128  points  with  a  time  bandwidth  product  NW  of  4 


Mean  and  envelope  of  the  64  spectral  densities 
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7.2.  MW  Estimate  On  128  points  With  NW  =  2 

512  spectra  were  computed  based  on  the  first  128  points  of  the  FFTs  576 
point  long  sequence.  The  time  bandwidth  product  NW  is  2,  the  number  of  win¬ 
dows  4  and  no  reshaping  is  performed.  The  Fourier  transforms  are  done  on  512 
points.  The  512  MW  spectra  are  overlaid  on  Figure  7.2.1  and  the  mean  and 
envelope  of  the  512  estimates  are  plotted  of  Figure  7.2.2. 

Because  the  averaging  MW  filter  normalized  bandwidth  is  only  0.03125, 
the  noise  dip  around  the  normalized  frequency  0.0625  (or  8  Hz)  is  resolved  and 
the  overall  shape  of  the  mean  spectrum  fits  the  theoretical  one.  Nevertheless, 
there  is  some  smearing  effect  at  DC  due  to  the  MW  averaging  filter  and  there  is 
much  more  scatter  than  for  the  FFT  estimates. 

The  random  error  e,  is  plotted  on  Figure  7.2.3.  It  is  higher  than  for  the 
periodogram  estimates  since  it  varies  around  0.6  at  high  frequency,  the  random 
error  has  a  maximum  of  16.9.  The  bias  error  tb  is  plotted  on  Figure  7.2.4  and  it 
is  also  slightly  higher  than  the  FFT  one,  although  bounded  by  0.2  at  high  fre¬ 
quencies  with  a  maximum  at  low  frequency  equal  to  13.26.  The  rms  error  £„*,  is 
therefore  higher  than  the  periodogram  rms  error  and  oscillates  around  0.6  at  high 
frequencies  as  shown  on  Figure  7.2.5.  The  rms  error  maximum  is  21.52  at  low 
frequency. 

Distributions  and  histograms  at  the  pre-selected  frequencies  (§2)  are  plotted 
on  Figures  7.2.6,  7.2.7,  7.2.8.  Since  the  number  of  windows  2NW  is  4,  one  can 
expect  a  number  of  degrees  of  freedom  smaller  or  equal  to  8.  The  distribution  at 
low  frequency  (10/512)  does  not  follow  a  chi-square  distribution  of  any  order. 
The  x22  surimposed  on  the  histogram  is  rejected  by  the  chi-square  goodness  of  fit 
test  at  the  a  =  0.05  level  of  significance.  The  MW  method  at  low  frequency  has  a 
poor  distribution.  One  can  estimate  the  number  of  degrees  of  freedom  to  be 
smaller  than  2,  which  means  that,  as  a  result  of  the  adaptive  weighting,  only  one 
eigenspectrum  is  used  in  the  MW  estimate. 

The  distribution  at  intermediate  frequency  (50/512)  follows  a  chi-square 
distribution  with  8  degrees  of  freedom  while  at  higher  frequency  (200/512)  the 
number  of  degrees  of  freedom  is  only  6. 
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Figure  7.2.1 


512  overlaid  spectral  densities  computed  by  the  MW  method  on 
128  points  and  a  time  bandwidth  product  NW  of  2 


Figure  7.2.2  Mean  and  envelope  of  the  512  spectral  densities 


Figure  7.2.3  Normalized  random  error  er 


Figure  7.2.4  Normalized  bias  error  £b 


Figure  7.2.5  Normalized  rms  error  e„ 


Distribution  af  the  Spectral  Estimate 
based  an  512  spectra _ _  4  /  4  /  86 


Figure  7.2.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  x22  distributions 


Figure  7.2.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x2»  distributions 


Figure  7.2.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  distributions 


7.3.  MW  Estimate  On  128  points  With  NW  =  3 

512  spectral  densities  were  computed,  as  in  §7.2  and  on  the  same  data,  with 
a  time  bandwidth  product  NW  of  3  and  six  windows.  No  reshaping  was  per¬ 
formed. 

The  overlaid  spectra  are  plotted  on  Figure  7.3.1.  The  mean  and  envelope  of 
the  estimates  are  plotted  on  Figure  7.3.2  Since  the  averaging  filter  normalized 
bandwidth  is  0.04687,  the  estimates  are  expected  to  have  more  bias  but  less  vari¬ 
ance  [2,3].  There  is  indeed  more  smearing  at  DC  than  before  and  the  shipping 
bump  or  power  spectrum  local  maximum  is  slightly  broader.  On  the  other  hand, 
one  can  easily  visualize  the  decrease  in  variance  with  respect  to  the  results  of 
§7.2,  where  NW  equal  2. 

The  random  error  er,  plotted  on  Figure  7.3.3,  is  essentially  varying  around 
0.5,  thus  lower  than  for  a  time  bandwidth  product  NW  of  2.  It  has  a  maximum 
value  of  26.  The  bias  error  ib,  plotted  on  Figure  7.3.4,  is  at  low  frequency  higher 
than  for  a  time  bandwidth  product  NW  of  2,  but  at  high  frequency  it  is  similar  to 
the  results  of  §6  (FFT  estimates).  The  maximum  value  of  th  at  low  frequency  is 
21.2.  The  rms  error,  plotted  on  Figure  7.3.5,  is  lower  than  for  a  time  bandwidth 
product  of  2  at  high  frequencies  :  0.5  instead  of  0.6.  At  low  frequencies,  where 
the  spectrum  is  red,  er  goes  up  to  33.55. 

The  distributions  of  the  spectral  estimates  at  the  selected  frequencies  (§2) 
are  plotted  on  Figure  7.3.6,  7.3.7  and  7.3.8.  The  number  of  degrees  of  freedom  is 
expected  to  be  smaller  or  equal  to  12  since  the  MW  estimate  is  an  average  of  six 
eigenspectra.  Like  in  §7.2,  the  low  frequency  (10/512)  histogram  does  not  follow 
a  chi-square  law  of  any  order.  The  x22  surimposed  on  the  histogram  (Figure  7.3.6) 
was  rejected  by  the  chi-square  goodness  of  fit  at  the  a  =  0.05  level  of  significance. 
Once  again,  the  distribution  is  poor  and  the  corresponding  number  of  degrees  of 
freedom  of  estimate  is  smaller  than  2. 

At  intermediate  and  high  frequency  (50/512  and  200/512),  the  sample  dis¬ 
tribution  is  fitted  by  a  x2w  ■ 
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Figure  7.3.3  Normalized  random  error  er 


Figure  7.3.4  Normalized  bias  error  eb 


Figure  7.3.5  Normalized  rms  error  e„ 
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Figure  7.3.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  x\  distributions 


Figure  7.3.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x\o  distributions 


Figure  7.3.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2io  distributions 


7.4.  MW  Estimate  On  576  Points  with  NW=4 

To  have  a  more  complete  picture  of  the  MW  estimation  capabilities,  512 
MW  estimates  based  on  the  same  576  points  as  in  §6  (JFT  estimates)  are  com¬ 
puted.  The  Fourier  transform  are  performed  on  2048  points. 

The  overlaid  spectra  are  plotted  on  Figure  7.4.1  and  the  mean  and  envelope 
of  the  estimate  are  plotted  on  Figure  7.4.2.  One  can  already  see  that  the  MW 
method,  in  this  case,  resolves  well  the  spectrum  with  a  variance  control  similar  to 
the  FFT  estimates  of  §6. 

The  random  error  er,  plotted  on  Figure  7.4.3,  is  essentially  equal  to  0.4  with 
a  peak  to  2.8  at  low  frequency.  The  bias  error  eh,  plotted  on  Figure  7.4.4,  varies 
between  0  and  0.1  with  an  extremum  down  to  -0.1  corresponding  to  the  local 
maximum  of  the  spectrum  (shipping  bump)  and  with  a  maximum  of  value  2.85  at 
low  frequency.  The  rms  error  e^,  is  dominated  by  the  random  error  as  can  be 
seen  on  Figure  7.4.5  ;  it  varies  around  0.4  with  a  peak  value  at  low  frequency  of 
4.  These  errors  are  comparable  to  the  errors  obtained  with  the  FFT  based  esti¬ 
mates  (§6). 

The  distributions  of  the  spectral  estimates  at  the  pre-selected  frequencies 
(§2)  are  plotted  on  Figure  7.4.6,  7.4.7,  7.4.8.  The  maximum  number  of  degrees 
of  freedom  one  can  expect  is  16,  since  8  windows  are  used  (NW=4).  As  before 
chi-square  distributions  are  fitted  to  the  sample  histograms  at  the  a =0.05  level  of 
significance,  the  numbers  of  degrees  of  freedom  are  at  the  normalized  frequen¬ 
cies  10/512,  50/512  and  200/512  respectively  12, 16  and  14. 
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Figure  7.4.1  512  overlaid  spectral  densities  computed  by  the  MW  method  on 

576  points  with  a  time  bandwidth  product  NW  of  4 


Figure  7.4.2 


Mean  and  envelope  of  the  spectral  densities 


Figure  7.4.3  Normalized  Random  error  e. 


Figure  7.4.4  Normalized  bias  error  e* 


Figure  7.4.5  Normalized  rms  error  e™ 
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Figure  7.4.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  %2n  distributions 


Figure  7.4.7 


Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  %2ie  distributions 


Figure  7.4.8 


Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2i4  distributions 
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8.  Maximum  Likelihood  Estimate 


8.1.  Order  Selection 

As  noted  in  §4,  the  Maximum  Likelihood  spectral  estimate  of  order  p  is  an 
average  of  the  AR  spectral  estimates  of  order  smaller  or  equal  to  p.  The  order 
selection  for  the  MLM  method  can  be  thought  similar  to  the  order  selection  of  an 
AR  spectral  estimate. 

A  way  to  select  an  AR  model  order  is  to  look  at  the  residual  energy  between 
the  data  and  the  modelled  data.  Such  residual  decreases  with  the  model  order  and 
eventually  settles.  Such  an  analysis  was  carried  out  with  128  data  points 
extracted  from  the  "ambient  noise"  data  set  :  the  residual  energy  is  plotted  on 
Figure  8.1.1. 

The  residual  seems  to  level  off  for  an  AR  model  order  of  16  which  should 
be  a  good  choice  for  the  maximum  likelihood  spectral  estimate.  To  have  more 
information  in  the  comparison  between  the  MW  estimates  and  the  MLM  esti¬ 
mates,  the  simulation  is  carried  out  for  orders  p  equal  to  8, 16  and  32. 
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8.2.  MLM  Estimate  of  order  8  On  128  points 

512  spectra  are  computed  on  the  same  128  data  points  as  the  MW  spectral 
estimates  of  §2  and  §3.  The  order  p  of  the  autocorrelation  matrix  Rp  is  8.  The 
data  are  well  behaved  since  none  of  the  autocorrelation  matrix  was  singular.  The 
Fourier  Transforms  that  implements  the  method  are  done  on  256  points. 

The  overlaid  spectral  density  realizations  obtained  by  the  MLM  method  are 
plotted  on  Figure  8.2.1.  The  mean  and  envelope  of  the  estimates  are  plotted  on 
Figure  8.2.2.  The  selection  of  a  low  order  produces  a  good  variance  control  (the 
spectral  estimates  are  very  smooth),  but  a  poor  resolution.  The  dip  at  the  normal¬ 
ized  frequency  0.0625  is  not  resolved  and  the  shape  of  the  spectrum  does  not  fit 
the  theoretical  one.  This  lack  of  resolution  could  be  expected  since  the  dip 
corresponds  to  the  moving  average  pan  of  the  ARMA  (5,4)  process  used  to 
create  the  "ambient  noise"  data.  There  is  a  large  scatter  of  the  estimates  at  low 
frequency. 

The  random  error  er,  plotted  on  Figure  8.2.3,  is  varying  below  0.3  at  inter¬ 
mediate  and  high  frequency  with  a  maximum  value  of  7.5  at  low  frequency.  The 
bias  error  eb,  plotted  on  Figure  8.2.4,  is  generally  below  0.2  except  at  low  fre¬ 
quency  where  it  goes  up  to  11.  The  rms  error  e,*,,  plotted  on  Figure  8.2.5,  stays 
below  0.3  in  most  the  frequency  band  with  a  peak  value,  at  the  power  spectrum 
local  maximum  (shipping  bump),  of  0.7.  The  highest  value  of  e**,  is  13.6  at  low 
frequency. 

The  distributions  and  sample  histograms  at  the  pre-selected  frequencies  (§2) 
are  plotted  on  Figures  8.2.6,  8.2.7  and  8.2.8.  The  low  frequency  (10/512)  histo¬ 
gram  look  like  the  ones  of  §7.2  or  §7.3  (MW  estimate  with  NW  equal  2  or  3). 
Any  chi-square  distribution  of  any  order  is  rejected  at  the  a  -0.05  level  of 
significance.  At  the  frequency  50/512,  the  fitted  chi-square  accepted  by  the 
goodness  of  fit  test  has  16  degrees  of  freedom.  At  high  frequency  (200/512),  the 
distribution  is  essentially  a  Gaussian  distribution  :  the  fitted  chi-square  has  36 
degrees  of  freedom. 


_ Normalize 

Characteristic  Value 
Maximum  Value 


Approximate  Number  of  Dee 
Normalized  Frequency  10/512 

Degrees  of  Freedom  <2 


rees  of  Freedom 
50/512  1  200/512 
16  36 


'  VaVAC'.V-V,V  ■  v 


■  %  v  * 


If)  o 


..ivvyj 


Figure  8.2.1 


512  overlaid  spectral  densities  computed  by  the  MLM  method 
on  128  points  with  an  order  p  of  8 


Figure  8.2.2  Mean  and  envelope  of  the  spectral  densities 


Figure  8.2.3  Normalized  random  error  tr 


Figure  8.2.4  Normalized  bias  error  £4 


Figure  8.2.5  Normalized  rms  error  e„ 
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Figure  8.2.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  distributions 


Figure  8.2.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x216  distributions 


Figure  8.2.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  y2m  distributions 
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8.3.  MLM  Estimate  of  order  16  On  128  Points 

The  order  of  the  autocorrelation  matrix  is  doubled,  p  is  16,  one  expects 
more  resolution  and  also  more  variance  than  in  §8.2.  512  spectral  estimates 
based  on  the  usual  128  points  of  the  "ambient  noise”  sequence  are  computed. 
The  Fourier  transforms  implementing  the  MLM  method  are  performed  on  256 
points. 

The  overlaid  spectra  are  plotted  on  Figure  8.3.1.  The  mean  and  envelope  of 
the  estimates  are  plotted  on  Figure  8.3.2.  One  can  already  recognize  an  improve¬ 
ment  in  the  resolution  of  the  estimate  with  respect  to  §8.2,  for  a  limited  variance 
increase.  There  is  still  a  large  scatter  of  the  spectral  estimates  near  DC. 

The  random  error  er,  plotted  on  Figure  8.3.3,  is  essentially  equal  to  0.3 
except  at  high  frequency  where  it  increases  up  to  0.4,  and  at  low  frequency 
where  it  peaks  up  to  3.5.  The  random  error  has  a  minimum  at  the  power  spec¬ 
trum  local  maximum  (shipping  bump)  with  a  value  of  0.2. 

The  bias  error  et,  plotted  on  Figure  8.3.4,  is  essentially  zero  with  two 
extrema,  the  first  one  at  low  frequency  with  value  equal  to  6.2,  and  the  second 
one  at  the  power  spectrum  local  maximum  (shipping  bump)  with  a  value  equal  to 
-0.6.  The  rms  error  is  plotted  on  Figure  8.3.5,  it  varies  as  in  §8.2  around  0.3 
except  at  high  frequency  where  it  increases  to  0.4  and  at  low  frequency  where  it 
has  a  maximum  value  of  7. 

The  distributions  and  sample  histograms  at  the  pre-selected  frequencies  (§2) 
are  plotted  on  Figure  8.3.6,  8.3.7,  8.3.8.  Chi-square  distributions  are  fitted  to  the 
sample  histograms  at  the  a  =  0.05  level  of  significance  with  a  number  of  degrees 
of  freedom  equal  to  6,  12  and  20  for  the  respective  normalized  frequencies 
10/512,  50/512  and  200/512.  Although  the  variance  only  increases  slightly,  the 
order  selection  has  an  impact  on  the  distributions. 


Table  8A1 

Normalized  Errors 

tr 

e„ 

Characteristic  Value 

0.3 

=  0 

0.3 

Maximum  Value 

3.5 

6.2 

7 

_ Table  8.3.2 

ADDroximate  Number  of  Deerees  of  Freedom 

Normalized  Freauencv 

10/512 

50/512 

200/512 

Decrees  of  Freedom 

6 

12 

20 

<& 


ft 


► 


5 

$ 

6 


I 

* 


7A 
$ 


.  A  » .  -va-. 


n 


A.  ■ 


Figure  8.3.1  512  overlaid  spectral  densities  computed  by  the  MLM  method 

on  128  points  with  an  order  p  equal  16 


Figure  8.3.2  Mean  and  envelope  of  the  spectral  densities 
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Figure  8.3.3  Normalized  random  error  e, 


Figure  8.3.4  Normalized  bias  error  e* 


Figure  8.3.5  Normalized  rms  error  e„ 
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Figure  8.3.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  %2s  distributions 


Figure  8.3.7 


Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  %\2  distributions 


Figure  8.3.8 


Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2»  distributions 
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8.4.  MLM  Estimate  Of  Order  32  On  128  Points 


The  order  p  of  the  autocorrelation  matrix  is  now  32.  Once  again,  the  512 
spectral  estimates  based  on  the  same  128  data  points  are  computed.  One  expects 
more  resolution  and  more  variance  than  in  §  8.3. 

The  overlaid  spectra  are  plotted  on  Figure  8.4.1.  The  mean  and  envelope  of 
the  estimates  are  plotted  on  Figure  8.4.2.  One  visually  checks  an  increase  in  the 
resolution  of  the  estimate  and  its  variability. 

The  random  error  ef,  plotted  on  Figure  8.4.3,  varies  below  0.4  except  at 
high  frequencies  where  it  increases  up  to  0.48  and  at  low  frequency  where  it 
peaks  up  to  1.6.  As  before  there  is  a  minimum  value  of  0.25  that  corresponds  to 
the  power  spectrum  local  maximum  (shipping  bump).  The  bias  error  tb ,  plotted 
on  Figure  8.4.4,  varies  around  -0.2  with  two  peaks,  one  of  value  2  at  low  fre¬ 
quency  and  one  of  value  -0.55  at  the  power  spectrum  local  maximum.  The  rms 
error  is  plotted  on  Figure  8.4.5,  it  has  a  peak  at  low  frequency  of  2.5  and  then 
varies  between  0.4  and  0.5. 

The  distributions  and  sample  histograms  are  plotted  on  Figures  8.4.6,  8.4.7 
and  8.4.8.  Although  the  increase  of  the  order  p  has  a  weak  impact  on  the  random 
error,  it  significantly  decreases  the  number  of  degrees  of  freedom  of  the  estimate. 
As  before  chi-square  distributions  were  fitted  using  the  chi-square  goodness  of  fit 
test  with  a  level  of  significance  a  =  0.05.  The  number  of  degrees  of  freedom  are 
respectively  6, 6, 10  for  the  normalized  frequencies  10/512, 50/512  and  200/512. 


Table  8.4.1 

Normalized  Errors 

e. 

e* 

e_. 

Characteristic  Value 

0.4 

-0.2 

<0.5 

Maximum  Value 

1.6 

2 

2.5 

_ Table  8A2 

Amjroximate  Number  of  Deerees  of  Freedom 

Normalized  Freauencv 

10/512 

50/512 

200/512 

Deerees  of  Freedom 

6 

6 

10 

Figure  8.4.1  512  overlaid  spectral  densities  computed  by  the  MLM  method 

on  128  points  with  an  order  p  of  32 


Figure  8.4.2  Mean  and  envelope  of  the  spectral  densities 


Figure  8.4.3  Normalized  random  error  e 


Figure  8.4.4  Normalized  bias  enor  e* 


Figure  8.4.5  Normalized  rms  error  e, 


Figure  8.4.6  Sample  histogTam  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  x26  distributions 


Figure  8.4.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x\  distributions 


Figure  8.4.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2io  distributions 


40 


8.5.  MLM  Estimate  Of  Order  32  On  576  points 

The  MLM  spectral  estimates  are  computed  on  the  same  576  points  as  in  §6 
and  §7.4.  More  data  points  are  used  in  the  estimation  of  the  autocorrelation 
matrix  of  order  32.  One  expects  more  resolution  and  better  variance  control. 

The  overlaid  spectra  are  plotted  on  Figure  8.5.1.  The  mean  and  envelope  of 
the  estimates  are  plotted  on  Figure  8.5.2.  One  observes  a  significant  decrease  in 
the  variability  of  the  estimate. 

The  random  error  er,  plotted  on  Figure  8.5.3  is  essentially  equal  to  0.2  with 
a  low  maximum  value  at  0.54.  The  bias  error  efc,  plotted  on  Figure  8.5.4,  is  essen¬ 
tially  0  with  two  extrema,  one  at  the  low  frequency  dip  and  of  value  1.9,  the 
other  at  the  power  spectrum  local  maximum  of  value  -0.4.  The  rms  error  e^,  is 
plotted  on  figure  8.5.5  and  varies  around  0.2  except  at  low  frequency  where  is 
has  a  peak  to  1.8,  with  a  secondary  maximum  of  value  0.4  at  the  power  spectrum 
local  maximum. 

The  distributions  and  sample  histograms  at  the  pre-selected  frequencies  (§2) 
are  plotted  on  Figures  8.5.6,  8.5.7,  8.5.8.  Comparison  with  the  fitted  Gaussian 
law  indicates  that  the  estimate  has  a  large  number  of  degrees  of  freedom.  At  low 
and  high  frequencies,  an  attempt  to  fit  a  chi-square  distribution  shows  that  the 
number  of  degrees  of  freedom  is  greater  than  40.  At  intermediate  frequency 
(50/512),  the  test  fit  is  achieved  for  a  x2x  (36  degrees  of  freedom).  The  goodness 
of  fit  test  is  nevertheless  rejected  at  the  a =0.05  level  of  significance.  The  particu¬ 
lar  statistical  realization  has  an  unlikely  bimodal  distribution  that  affects  the  chi- 
square  goodness  of  fit  test  at  this  particular  frequency. 


Figure  8.5.1  512  overlaid  spectral  densities  computed  by  the  MLM  method 

on  576  points  with  an  order  p  of  32 


Figure  8.5.2  Mean  and  envelope  of  the  spectral  densities 


Figure  8.5.3  Normalized  random  error 


Figure  8.5.4  Normalized  bias  error  tb 


Figure  8.5.5  Normalized  rms  error  e_. 


1% 
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Figure  8.5.6  Sample  histogram  at  the  normalized  frequency  10/512  overlaid 
with  fitted  Gaussian  and  x24o  distributions. 


Figure  8.5.7  Sample  histogram  at  the  normalized  frequency  50/512  overlaid 
with  fitted  Gaussian  and  x2%  distributions. 


Figure  8.5.8  Sample  histogram  at  the  normalized  frequency  200/512  over¬ 
laid  with  fitted  Gaussian  and  x2«  distributions 
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9.  Conclusions 

The  results  are  summarized  in  terms  of  the  normalized  error  terms  (§2)  in 
Table  9.1,  9.2,  9.3  and  in  terms  of  the  number  of  degrees  of  freedom  in  Table 
9.4.  The  reference  estimator  is  the  conventional  FFT  method  based  on  576 
points  using  8  incoherently  averaged  FFTs  on  128  data  points  with  50  %  overlap 
and  a  Hanning  window.  The  tables  include  a  relative  efficiency  measure  defined 
as  the  percentage  ratio  of  the  rms  or  random  error  of  the  reference  estimator 
(averaged  FFTs)  to  the  rms  or  random  error  of  the  particular  method  [9]. 

The  MLM  spectral  estimates  of  order  32  and  the  MW  spectral  estimates 
with  a  time  bandwidth  NW  of  4,  based  on  the  same  576  data  sequences,  allow  a 
comparison  of  the  three  methods  for  a  fixed  number  of  points.  In  all  cases,  the 
underlying  spectral  density  is  well  resolved.  The  MLM  estimator  seems  to  be  the 
best  in  this  simulation  when  bias,  random  and  rms  errors  are  considered. 

The  resolution  of  the  MW  method  when  processing  576  points  is  given  by 
the  Slepian  window  bandwidth  2W,  in  the  present  case  2W  is  8/576  or  1/72.  The 
resolution  of  the  conventional  FFT  based  method  is  given  by  the  Hanning  win¬ 
dow  width  which  is  approximately  two  bins  [10]  or  2/128  that  is  1/64.  The 
results  given  by  these  two  methods  are  equivalent  in  term  of  the  normalized 
errors,  with  a  barely  better  conventional  estimate. 

The  distributions  of  the  MLM  estimates  of  order  32  and  on  576  points  are 
almost  Gaussian.  On  the  contrary,  the  MW  and  FFT  based  distributions  are 
almost  equivalent  chi-squares,  although  the  MW  method  seems  to  lose  four 
degrees  of  freedom  at  low  frequency  (in  the  red  and  low  magnitude  part  of  the 
power  spectral  density)  with  respect  to  the  conventional  FFT  estimates.  This  is 
due  to  the  adaptive  weighting  of  the  eigenspectra. 

Varying  statistical  properties  with  the  spectral  density  shape  and  dynamic 
range  are  also  characteristic  of  the  MLM  estimators  when  the  number  of  points 
in  the  data  sequence  is  reduced  to  128.  Such  variations  are  dramatic,  given  the 
more  stringent  conditions.  The  MW  method  and  the  MLM  method  of  order  8  on 
128  data  points  loose  statistical  consistency  at  low  frequency  in  the  neighbor¬ 
hood  of  the  power  spectral  density  dip  near  the  normalized  frequency  0.0625. 

In  this  simulation  on  short  data  sequence,  the  results  seem  also  favorable  to 
the  MLM  method  of  order  16  or  32.  In  addition  to  good  statistical  consistency, 
that  is  a  larger  number  of  degrees  of  freedom,  especially  at  low  frequency,  the 
average  rms  and  random  error  are  lower  than  the  MW  method  ones.  This  is  espe¬ 
cially  clear  looking  at  the  relative  efficiencies  with  respect  to  the  conventional 
FFT  method. 

At  high  frequency,  the  rms  and  random  error  efficiency  of  the  MLM 
method,  for  128  point  long  sequence  and  an  order  equal  to  16  or  32,  is  above 
100%  ;  this  is  also  the  case  when  the  MLM  estimates  are  computed  on  576 
points.  It  means  that  the  estimates  are  better  than  the  conventional  FFT  estima¬ 
tor.  This  result  is  expected  since  the  MLM  method  can  be  shown  to  be  equivalent 
to  an  average  of  AR  based  estimates.  Since  data  were  created  by  an  ARMA  pro¬ 
cess,  they  are  well  modeled  by  the  MLM  method  in  the  high  end  of  the  fre¬ 
quency  band.  The  general  results  are  not  likely  an  artifact  due  to  the  data  since 
the  low  frequency  end  (dip)  corresponds  to  the  moving  average  part  of  the  spec¬ 
trum.  Even  at  low  frequency,  the  errors  are  lower  and  the  number  of  degrees  of 
freedom  is  higher  in  the  MLM  estimates  with  p  equal  16  or  32  than  for  the  MW 
estimates. 
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